#delimit;
clear;

use 20161021_Working_hh.dta, clear;

replace Y=(Y/1000)*CPIU2016;
replace totalTransfers=(totalTransfers/1000)*CPIU2016;

statsby _b _se n=(e(N)), by(year) verbose nodots clear: 
 nl (totalTransfers = {a} + {b1}*exp({b2}*Y)), initial(a 1 b1 2 b2 -.1) vce(gnr);

rename a_b A;
rename b1_b B1;
rename b2_b B2; 


gen cov=0;
replace a_se=a_se^2;
replace b1_se=b1_se^2;
replace b2_se=b2_se^2;

rename a_se V_A;
rename b1_se V_B1;
rename b2_se V_B2;

rename _eq7_n n;

save 20161021_parametersByYear, replace;



*>>>>A  Constructing Matrices;

use 20161021_parametersByYear, clear;

matrix drop _all;
sort year;
mkmat A B1 B2, matrix(BETA);
matrix rownames BETA=  1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
       2008 2009 2010 2011 2012 2013 2014 2015
; 


mkmat V_A cov V_B1 cov cov V_B2  , matrix(COV);
matrix rownames COV= 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
       2008 2009 2010 2011 2012 2013 2014 2015;


local list "1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007
	          2008 2009 2010 2011 2012 2013 2014 2015"; 


foreach i of local list  {;

    matrix BETA_`i'= BETA["`i'", 1..3];
    matrix COV_`i'=COV["`i'", 1..6];
        
};


*>>>>B  Generating Simulated Parameters;

drop _all;
save "20161021_SimulationsAnnual.dta", replace emptyok;

foreach i of local list  {;


    drawnorm ALPHA BETA_1 BETA_2, n(1000) means(BETA_`i') cov(COV_`i') cstorage(lower) seed(00001) clear;
    generate year=`i';
    
	append using "20161021_SimulationsAnnual.dta";
    save "20161021_SimulationsAnnual.dta", replace;
 
    
};
*/;
*>>>>C  Generating Distribution of TAU and R;


use 20161021_SimulationsAnnual.dta, clear;
sort year;

merge m:1 year using 20161021_parametersByYear; 



gen PSI_1=.5*19.790;
gen PSI_2=1*19.790;
gen PSI_3=1.5*19.790;
gen PSI_4=2*19.790;



forvalues p=1/4{;

generate TAU_`p'=.7*PSI_`p';

foreach i of local list  {;

local j=1;

while `j'  < 100  {;
    generate TAU_`p'_`i'_`j'= TAU_`p' +`j'*.1 if (year==`i') ;
	generate TEST_`p'_`i'_`j'=(PSI_`p'-TAU_`p'_`i'_`j') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`j')) if (year==`i') ;
	replace TAU_`p'=TAU_`p'_`i'_`j' if TEST_`p'_`i'_`j'<.1 & TEST_`p'_`i'_`j'>0;
	
	local j=`j'+1;
	};



local k=1;

while `k'  < 100  {;
    replace TAU_`p'_`i'_`k'= TAU_`p' - `k'*.01 if (year==`i') ;
	replace TEST_`p'_`i'_`k'=(PSI_`p'-TAU_`p'_`i'_`k') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`k')) if (year==`i') ;
	replace TAU_`p'=TAU_`p'_`i'_`k' if TEST_`p'_`i'_`k'<.01 & TEST_`p'_`i'_`k'>0;
	
	local k=`k'+1;
	};


local l=1;

while `l'  < 100  {;
    replace TAU_`p'_`i'_`l'= TAU_`p' - `l'*.001 if (year==`i') ;
	replace TEST_`p'_`i'_`l'=(PSI_`p'-TAU_`p'_`i'_`l') - (ALPHA+BETA_1*exp(BETA_2*TAU_`p'_`i'_`l')) if (year==`i') ;
	replace TAU_`p'=TAU_`p'_`i'_`l' if TEST_`p'_`i'_`l'<.001 & TEST_`p'_`i'_`l'>0;
	
	local l=`l'+1;
	};



drop TEST* TAU_`p'_*;

};

generate R_`p'= ((A*TAU_`p' + (B1/B2)*exp(B2*TAU_`p') - (B1/B2))  + (PSI_`p'^2 - PSI_`p'^2/2) - (PSI_`p'*TAU_`p' - TAU_`p'^2/2)) /(PSI_`p'^2 - PSI_`p'^2/2);


};
                                                                                  

/* NB Corrected equation */;     


#delimit;
collapse A B1 B2  V_A V_B1 V_B2 TAU* PSI* R* n (sd) SD_R_1=R_1 SD_R_2=R_2 SD_R_3=R_3 SD_R_4=R_4 SD_TAU_2=TAU_2, by(year); 


save "20161021_Annual_Parameters_w_Var_8814.dta", replace;






